Analysis of separate coil data (chris) white matter

In [1]: 
import numpy as np
import numpy.linalg as nla
import matplotlib.pyplot as plt
import nibabel as nib
import os
%matplotlib inline
In [2]: 
path = '/home/snarles/predator'
In [4]: 
os.listdir(path)
#os.listdir(path + '/8631_5_1_pfile/coil_images')
Out[4]:
['8631_5_1_pfile',
 '8631_2_1_wm_mask.nii.gz',
 'chris2_bvec.csv',
 'chris1_bvec.csv',
 '8631_11_1_pfile']
In [6]: 
bvecs = np.loadtxt(path + '/chris1_bvec.csv', delimiter = ' ').T
b0_inds = list(np.where(np.array([nla.norm(v) for v in bvecs]) == 0)[0])
b0_inds = b0_inds[2:]
#plt.scatter(np.array(b0_inds),np.ones(10))
In [76]: 
#wm = nib.load(path + '/8631_2_1_wm_mask.nii.gz').get_data()
wm = np.ones((120, 120, 69))
In [9]: 
np.shape(wm)
Out[9]:
(120, 120, 69)
In [77]: 
c1_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil18_ec.nii.gz').get_data()
c1_0 = c1_0[:, :, :, b0_inds]
c2_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil21_ec.nii.gz').get_data()
c2_0 = c2_0[:, :, :, b0_inds]
c3_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil28_ec.nii.gz').get_data()
c3_0 = c3_0[:, :, :, b0_inds]
c4_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil31_ec.nii.gz').get_data()
c4_0 = c4_0[:, :, :, b0_inds]
In [78]: 
c1 = c1_0
c2 = c2_0
c3 = c3_0
c4 = c4_0
for i in range(10):
    c1[:, :, :, i] = c1_0[:, :, :, i] * wm
    c2[:, :, :, i] = c2_0[:, :, :, i] * wm
    c3[:, :, :, i] = c3_0[:, :, :, i] * wm
    c4[:, :, :, i] = c4_0[:, :, :, i] * wm
c1_mu = np.mean(c1, axis = 3)
c2_mu = np.mean(c2, axis = 3)
c3_mu = np.mean(c3, axis = 3)
c4_mu = np.mean(c4, axis = 3)
c1_0mu = np.mean(c1_0, axis = 3)
c2_0mu = np.mean(c2_0, axis = 3)
c3_0mu = np.mean(c3_0, axis = 3)
c4_0mu = np.mean(c4_0, axis = 3)
In [25]: 
(np.min(c1_mu[:]), np.max(c1_mu[:])), (np.min(c2_mu[:]), np.max(c2_mu[:])), (np.min(c3_mu[:]), np.max(c3_mu[:]))
Out[25]:
((-0.551247, 13.458018), (-0.22318967, 4.8017988), (-0.38862443, 4.839869))
In [26]: 
vmin_ = -1.0
vmax_ = 15.0
figsize_ = (3, 3)
z_inds = [10, 20, 30, 40, 50]

Sep coil images

In [79]: 
#w_inds = range(10)
#z_inds = [12, 14, 16, 18, 20]
figsize_2 = (len(z_inds) * figsize_[0], figsize_[1])
f_c1, axarr = plt.subplots(1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
    z_ind = z_inds[j]
    axarr[j].imshow(c1_mu[:, :, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
f_c2, axarr = plt.subplots(1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
    z_ind = z_inds[j]
    axarr[j].imshow(c2_mu[:, :, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
f_c3, axarr = plt.subplots(1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
    z_ind = z_inds[j]
    axarr[j].imshow(c3_mu[:, :, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
f_c4, axarr = plt.subplots(1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
    z_ind = z_inds[j]
    axarr[j].imshow(c4_mu[:, :, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)

shared variation

time is left-right

In [80]: 
w_inds = range(10)
vmin_2 = -1.0
vmax_2 = 1.0
figsize_2 = (len(w_inds) * figsize_[1]/2, 4 * figsize_[0])
f_sv, axarr = plt.subplots(4, len(w_inds)+1, sharex=True, sharey=True, figsize = figsize_2)
z_ind = 20
axarr[0, 0].imshow(c1_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[1, 0].imshow(c2_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[2, 0].imshow(c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[3, 0].imshow(c4_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
for i in range(len(w_inds)):
    w_ind = w_inds[i]
    axarr[0, i+1].imshow(c1[:, 20:80, z_ind, w_ind] - c1_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
    axarr[1, i+1].imshow(c2[:, 20:80, z_ind, w_ind] - c2_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
    axarr[2, i+1].imshow(c3[:, 20:80, z_ind, w_ind] - c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
    axarr[3, i+1].imshow(c4[:, 20:80, z_ind, w_ind] - c4_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)

z-slice effects

In [81]: 
z_inds = range(18, 26)
w_inds = range(10)
figsize_2 = (len(z_inds) * figsize_[0]/2, (len(w_inds)+1) * figsize_[1])
f_c1, axarr = plt.subplots(len(w_inds)+1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
    z_ind = z_inds[j]
    axarr[0, j].imshow(c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
for i in range(len(w_inds)):
    for j in range(len(z_inds)):
        w_ind = w_inds[i]
        z_ind = z_inds[j]
        axarr[i+1, j].imshow(c3[:, 20:80, z_ind, w_ind] - c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)

what motion artifacts would look like

In [82]: 
z_inds = range(18, 26)
figsize_2 = (len(z_inds) * figsize_[0]/2, 4 * figsize_[1])
f_c1, axarr = plt.subplots(4, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
vmin_3 = -2.0
vmax_3 = 2.0
for j in range(len(z_inds)):
    z_ind = z_inds[j]
    axarr[0, j].imshow((c3_0mu[10:100, 20:80, z_ind]) * wm[10:100, 20:80, z_ind],
                       cmap='bone', vmin = vmin_, vmax = vmax_)
    axarr[1, j].imshow((c3_0mu[10:100, 20:80, z_ind] - c3_0mu[9:99, 20:80, z_ind]) * wm[10:100, 20:80, z_ind],
                       cmap='bone', vmin = vmin_3, vmax = vmax_3)
    axarr[2, j].imshow((c3_0mu[10:100, 20:80, z_ind] - c3_0mu[10:100, 19:79, z_ind]) * wm[10:100, 20:80, z_ind],
                       cmap='bone', vmin = vmin_3, vmax = vmax_3)
    axarr[3, j].imshow((c3_0mu[10:100, 20:80, z_ind] - c3_0mu[10:100, 20:80, z_ind-1]) * wm[10:100, 20:80, z_ind],
                       cmap='bone', vmin = vmin_3, vmax = vmax_3)

Aggregating diffusion data

In [83]: 
b0_breaks = np.array([2] + b0_inds)
break_starts = b0_breaks[:-1]+1
break_ends = b0_breaks[1:]-1
break_lens = break_ends - break_starts
len(break_lens)
Out[83]:
10
In [84]: 
def diffusion_aggregate(data_0):
    data_d = np.zeros((120, 120, 69, 10))
    for i in range(10):
        data_d[:, :, :, i] = np.mean(data_0[:, :, :, break_starts[i]:break_ends[i]], axis = 3) * wm
    return data_d

def demean(data_0):
    data_0 = np.array(data_0)
    wl = np.shape(data_0)[3]
    data_mu = np.mean(data_0, axis = 3)
    for i in range(wl):
        data_0[:, :, :, i] = data_0[:, :, :, i] - data_mu
    return data_0

def da_dac(data_0):
    data_da = diffusion_aggregate(data_0)
    data_dac = demean(data_da)
    return data_da, data_dac
In [85]: 
c1_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil18_ec.nii.gz').get_data()
c1_da, c1_dac = da_dac(c1_0)
del c1_0
c2_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil21_ec.nii.gz').get_data()
c2_da, c2_dac = da_dac(c2_0)
del c2_0
c3_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil28_ec.nii.gz').get_data()
c3_da, c3_dac = da_dac(c3_0)
del c3_0
c4_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil31_ec.nii.gz').get_data()
c4_da, c4_dac = da_dac(c4_0)
del c4_0

shared variation

time is left-right

In [86]: 
w_inds = range(10)
vmin_2 = -1.0
vmax_2 = 1.0
figsize_2 = (len(w_inds) * figsize_[1]/2, 4 * figsize_[0])
f_sv, axarr = plt.subplots(4, len(w_inds)+1, sharex=True, sharey=True, figsize = figsize_2)
z_ind = 20
axarr[0, 0].imshow(c1_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[1, 0].imshow(c2_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[2, 0].imshow(c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[3, 0].imshow(c4_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
for i in range(len(w_inds)):
    w_ind = w_inds[i]
    axarr[0, i+1].imshow(c1_dac[:, 20:80, z_ind, w_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
    axarr[1, i+1].imshow(c2_dac[:, 20:80, z_ind, w_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
    axarr[2, i+1].imshow(c3_dac[:, 20:80, z_ind, w_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
    axarr[3, i+1].imshow(c4_dac[:, 20:80, z_ind, w_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
In []: